
	clear
	use Data_20170719.dta
	set seed 5500

		/* 
		Note that the marginal effect is set to increase from 500 to 1500 troops.
		This represents a change from 500 below the median non-zero value of UN deployments (~1000 troops) 
		to 500 above the median. 500 is also approximately the size of one US battalion. So we're increasing the deployment
		from 1 to 3 battalions.
		US Models are set to match this change.
		*/
		
		* Marginal Effect of an increase in peacekeepers across number of contributing countries
			levelsof contribcount, local(levels)	
			xtset ccode year
			probit pdest c.contribcount##c.ln_totalpeacekeepers intensity_high w_intensity_high w_psource w_pdest polity2 ln_pop ln_gdpwb  border, robust 
			
			sum ln_totalpeacekeepers if e(sample) & ln_totalpeacekeepers>0

			sum contribcount if e(sample)
			local min = r(min)
			local max = r(max)
			local interval = ((`max'-`min')/100)
			
			margins, at(ln_totalpeacekeepers=6.22 ln_totalpeacekeepers=7.32 contribcount=(`min'(1)`max') (mean) w_psource w_pdest polity2 ln_pop ln_gdpwb border w_intensity_high intensity_high=0) mcompare(bonferroni) vsquish level(90) post  		
			
			local a = 4
			local b = 3
			#delimit ;
			matrix list e(b);
			mlincom 2-1, clear stats(all)	;
			while `a'<=170 { ;
			mlincom `a'-`b', stats(all) add ;
			local a=`a'+2 ;
			local b = `b'+2 ;
			}
			;
			matrix list _mlincom			;
			svmat double _mlincom, names(col);
			range xvar `min' `max' 85 ;
			gen yline = 0 ;
			
			#delimit ;
			twoway rarea ll ul xvar,  color(gs12%50) lcolor(gs12%0) fi(100)
			|| line lincom xvar, lpattern(solid) msym(O) lwidth(.3) 
			yline(0)
			xlabel(#12, format(%9.0fc)) 
			ylabel(, format(%9.2fc))
			ytitle("Change in Predicted Probability")
			legend(off)
			xtitle("Troop Contributing Countries")
			subtitle("A")
			aspect(1)
			saving(me_pk_contribcount, replace) 
			;	

			#delimit cr
			
			graph export "me_pdest_contribcount.pdf", replace

				
		* Marginal Effect of an increase in peacekeepers across duration of pko presence.
		clear
		use Data_20170719.dta
		set seed 5500

			levelsof duration_pk_totalpeacekeepers, local(levels)	
			xtset ccode year
			probit pdest c.duration_pk_totalpeacekeepers##c.ln_totalpeacekeepers intensity_high w_intensity_high w_psource w_pdest   polity2 ln_pop ln_gdpwb border, robust 

			sum ln_totalpeacekeepers if e(sample) & ln_totalpeacekeepers>0
			
			sum duration_pk_totalpeacekeepers if e(sample)
			local min = r(min)
			local max = r(max)
			local interval = ((`max'-`min')/30)
			
			margins, at(ln_totalpeacekeepers=6.22 ln_totalpeacekeepers=7.32 duration_pk_totalpeacekeepers=(`min'(1)`max') (mean) w_psource w_pdest polity2 ln_pop ln_gdpwb border w_intensity_high intensity_high=0) mcompare(bonferroni) vsquish level(90) post  		
			
			local a = 4
			local b = 3
			#delimit ;
			matrix list e(b);
			mlincom 2-1, clear stats(all)	;
			while `a'<=44 { ;
			mlincom `a'-`b', stats(all) add ;
			local a=`a'+2 ;
			local b = `b'+2 ;
			}
			;
			matrix list _mlincom			;
			svmat double _mlincom, names(col);
			range xvar `min' `max' 15 ;
			gen yline = 0 ;
			
			#delimit ;
			twoway rarea ll ul xvar,  color(gs12%50) lcolor(gs12%0) fi(100)
			|| line lincom xvar, lpattern(solid) msym(O) lwidth(.3) 
			yline(0)
			xlabel(#12, format(%9.0fc)) 
			ylabel(, format(%9.2fc))
			xtitle("Duration of Peacekeeping Presence")
			ytitle("Change in Predicted Probability")
			legend(off)
			title("")
			subtitle("B")
			aspect(1)
			saving(me_pk_duration, replace)
			;
			#delimit cr
			
			graph export "me_pdest_duration.pdf", replace
		
			graph combine me_pk_contribcount.gph me_pk_duration.gph, cols(2)
			graph export "me_pk_conditional.pdf", replace
		


		* Marginal Effect of an increase in peacekeepers across number of contributing countries
		clear
		use Data_20170719.dta
		set seed 5500
			
			levelsof contribcount, local(levels)	
			xtset ccode year
			probit pdest c.contribcount##c.ln_totalpeacekeepers intensity_high w_intensity_high w_psource w_pdest polity2 ln_pop ln_gdpwb  border if totalpeacekeeper>0, robust 
			
			sum ln_totalpeacekeepers if e(sample) & ln_totalpeacekeepers>0
			
			sum contribcount if e(sample)
			local min = r(min)
			local max = r(max)
			local interval = ((`max'-`min')/100)
			margins, at(ln_totalpeacekeepers=6.22 ln_totalpeacekeepers=7.32 contribcount=(`min'(1)`max') (mean) w_psource w_pdest polity2 ln_pop ln_gdpwb border w_intensity_high intensity_high=0) mcompare(bonferroni) vsquish level(90) post  		
			
			local a = 4
			local b = 3
			#delimit ;
			matrix list e(b);
			mlincom 2-1, clear stats(all)	;
			while `a'<=166 { ;
			mlincom `a'-`b', stats(all) add ;
			local a=`a'+2 ;
			local b = `b'+2 ;
			}
			;
			matrix list _mlincom			;
			svmat double _mlincom, names(col);
			range xvar `min' `max' 85 ;
			gen yline = 0 ;
			
			#delimit ;
			twoway rarea ll ul xvar,  color(gs12%50) lcolor(gs12%0) fi(100)
			|| line lincom xvar, lpattern(solid) msym(O) lwidth(.3) 
			yline(0)
			xlabel(#12, format(%9.0fc)) 
			ylabel(, format(%9.2fc))
			ytitle("Change in Predicted Probability")
			legend(off)
			xtitle("Troop Contributing Countries")
			subtitle("C")
			aspect(1)
			saving(me_pk_contribcount_limited, replace) 
			;	
			#delimit cr
			
			graph export "me_pdest_contribcount_limited.pdf", replace
		

		* Marginal Effect of an increase in peacekeepers across duration of pko presence.
		clear
		use Data_20170719.dta
		set seed 5500

			levelsof duration_pk_totalpeacekeepers, local(levels)	
			xtset ccode year
			probit pdest c.duration_pk_totalpeacekeepers##c.ln_totalpeacekeepers intensity_high w_intensity_high w_psource w_pdest   polity2 ln_pop ln_gdpwb border if totalpeacekeeper>0, robust 

			sum ln_totalpeacekeepers if e(sample) & ln_totalpeacekeepers>0
			
			sum duration_pk_totalpeacekeepers if e(sample)
			local min = r(min)
			local max = r(max)
			local interval = ((`max'-`min')/30)
			
			margins, at(ln_totalpeacekeepers=6.22 ln_totalpeacekeepers=7.32 duration_pk_totalpeacekeepers=(`min'(1)`max') (mean) w_psource w_pdest polity2 ln_pop ln_gdpwb border w_intensity_high intensity_high=0) mcompare(bonferroni) vsquish level(90) post  		
			
			local a = 4
			local b = 3
			#delimit ;
			matrix list e(b);
			mlincom 2-1, clear stats(all)	;
			while `a'<=35 { ;
			mlincom `a'-`b', stats(all) add ;
			local a=`a'+2 ;
			local b = `b'+2 ;
			}
			;
			matrix list _mlincom			;
			svmat double _mlincom, names(col);
			range xvar `min' `max' 15 ;
			gen yline = 0 ;
			
			#delimit ;
			twoway rarea ll ul xvar,  color(gs12%50) lcolor(gs12%0) fi(100)
			|| line lincom xvar, lpattern(solid) msym(O) lwidth(.3) 
			yline(0)
			xlabel(#12, format(%9.0fc)) 
			ylabel(, format(%9.2fc))
			xtitle("Duration of Peacekeeping Presence")
			ytitle("Change in Predicted Probability")
			legend(off)
			title("")
			subtitle("D")
			aspect(1)
			saving(me_pk_duration_limited, replace)
			;
			#delimit cr
		
			graph combine me_pk_contribcount.gph me_pk_duration.gph me_pk_contribcount_limited.gph me_pk_duration_limited.gph, cols(2) imargin(1 1 1 1) xsize(12) ysize(10) scale(.75) ycommon
			
			graph export "me_pdest_conditional_fourpanel.pdf", replace

		
			
